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\ We investigate the adiabatic orbital evolution of a point particle in the Kerr spacetime 

due to the emission of gravitational waves. In the case that the timescale of the orbital 
• evolution is enough smaller than the typical timescale of orbits, the evolution of orbits is 

^ ' characterized by the change rates of three constants of motion, the energy E, the azimuthal 

, angular momentum L, and the Carter constant Q. For E and L, we can evaluate their change 

'/^ . rates from the fluxes of the energy and the angular momentum at infinity and on the event 

' horizon according to the balance argument. On the other hand, for the Carter constant, we 

' cannot use the balance argument because we do not know the conserved current associated 

, with it. Recently, Mino proposed a new method of evaluating the averaged change rate of the 

Carter constant by using the radiative field. In our previous paper we developed a simplified 
' scheme for practical evaluation of the evolution of the Carter constant based on the Mino's 

^ , proposal. In this paper we describe our scheme in more detail, and derive explicit analytic 

formulae for the change rates of the energy, the angular momentum and the Carter constant. 

U '. 

§1. Introduction 

> ■ 

^ . It is believed that supermassive black holes (SMBHs) reside in central nuclei of 

^ I many galaxies, and they occasionally capture a stellar mass compact object (SMCO) 

■ - - ' which surrounds them. Gravitational waves from such binary systems with extreme 

mass ratios bring us information on the orbits of SMCOs and the spacetime structure 
near black holes. Therefore such systems are considered to be one of the most 
important targets of the LISA space-based gravitational wave detector.^) In order 
to detect gravitational waves emitted by extreme mass ratio inspirals (EMRIs) and 
to extract physical information from them efficiently, we need to predict accurate 
theoretical waveforms in advance. Our goal along the line of this paper is to precisely 
calculate theoretical waveforms from EMRIs. 

To investigate gravitational waves from EMRIs, our strategy is to adopt the 
black hole perturbation method:^^ we consider metric perturbations induced by a 
SMCO in a black hole spacetime governed by a SMBH. We also assume that a 
SMCO is described by a point particle, neglecting its internal structure. Under 
the above approximations, we can calculate the metric perturbation evaluated at 
infinity to predict gravitational waveforms. At the lowest order with respect to the 
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mass ratio, we may calculate the metric perturbation by approximating the particle's 
orbit by a background geodesic. To step further, we consider the orbital shift from 
the background geodesic by taking account of the self-force induced by the particle 
itself. 

In a Schwarzschild background, we can assume that the orbit is in the equatorial 
plane from the symmetry without loss of generality. Hence, the orbital velocity 
can be specified solely by the energy and the azimuthal angular momentum of the 
particle. Namely, we can evaluate the orbital evolution from the change rates of the 
energy and the angular momentum. Their averaged change rates can be evaluated by 
using the balance argument; the energy and the angular momentum that a particle 
loses are equal to the ones that are radiated to the infinity or across the horizon 
as gravitational waves because of the conservation laws. In the limit of a large 
mass ratio, averaged change rates will be sufficient to determine the leading order 
effects on the orbital evolution due to the self-force. In this sense the leading order 
effects can be read from the asymptotic behavior of the metric perturbation in the 
Schwarzschild case. 

On the other hand, the third constant of motion, i.e., the Carter constant, is 
necessary in addition to the energy and the azimuthal angular momentum to specify 
a geodesic in a Kerr background. However, there is no known conserved current com- 
posed of gravitational waves that is associated with the Carter constant, and hence 
we cannot use the balance argument to evaluate the change rate of the Carter con- 
stant. Therefore we have to calculate the self-force acting on a particle. ^^'^^ When we 
calculate the self-force, we are faced with the regularization problem. Although the 
formal expression for the regularized self-force had been derived, ^^"''^ doing explicit 
calculation is not so straightforward. 

Gal'tsov®^ proposed a method of calculating the loss rates of the energy and 
angular momentum of a particle by using the radiative part of metric perturbation, 
which was introduced earlier by Dirac.^^ The radiative field is defined by half re- 
tarded field minus half advanced one, which is a homogeneous solution of the field 
equation. It was shown that the time-averaged loss rates of the energy and angu- 
lar momentum evaluated by using the radiative field are identical with the results 
obtained from the balance argument. Recently, Mino proved that the Gal'tsov's 
scheme also gives the correct averaged change rate of the Carter constant.^'') The 
Gal'tsov-Mino method has a great advantage that we do not need any regularization 
procedure because the radiative field is free from divergence from the beginning. 
In Ref. 11), we briefly reported that the formula for the adiabatic evolution of the 
Carter constant based on Gal'tsov-Mino method can be largely simplified. In this 
paper, we explain the derivation of this new formula in detail. Applying our new 
formula, we explicitly calculate the change rate of the Carter constant for orbits with 
small eccentricities and inclinations. 

This paper is organized as follows. In Sec. [21 we give a brief review of the 
Kerr geometry and the geodesic motion. Next, we show a practical prescription to 
calculate the time-averaged change rates of the constants of motion in Sec. 01 We 
also derive a simplified expression for the change rate of the Carter constant. In 
Sec. ^ applying our prescription, we calculate the change rates of the constants 
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of motion and then show the analytic formulae of them for slightly eccentric and 
inclined orbits. Finally we devote Sec. [3 to summarize this paper. In Appendix 
we show the derivation of the radiative part of metric perturbation. And we also 
give short reviews on analytic methods of solving the radial Teukolsky equation and 
obtaining the spheroidal harmonics in Appendices IbI and IHl 



§2. Geodesic motion in the Kerr spacetime 



In this section, we give a brief review on geodesies in the Kerr geometry. The 
metric of the Kerr spacetime in the Boyer-Lindquist coordinates is 

ds^ = - il — (ir dtdif + —dr^ 



U J S A 

+Ude^ + (^r^ + a^ + sin2 9^ sin^ Odip^ (2-1) 

where 

S = r'^ + a^cos^9, A = r^-2Mr + a^. 

M and aM are the mass and angular momentum of the black hole, respectively. 
There are two Killing vectors reflecting the stationary and axisymmetric properties 
of the Kerr geometry: 

Cf,) = (1,0,0,0), ef^) = (0,0,0,1). (2-2) 

In addition, the Kerr spacetime possesses a Killing tensor, 

K^,i, = 2Z'Z(^n^) + r'^g^^, (2-3) 

which satisfies Kf^^^.p-^ = 0, where the parenthese operating on the indices is the 
notation for symmetric part of tensors. Here we have introduced null vectors, 

m'' := ^ f msin6',0, 1,^— ) . (2-4) 

We consider a point particle moving in the Kerr geometry: 

z''[r) = {t,{r),r,[r),e,[r),ip,[r)), 
where r is the proper time along the orbit. Here we introduce quantities defined by 

E := = (l - ^) «* + ^J^^^^^:f^u'. (2.5) 

sm Uz 
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where = dz^/dr. These quantities remain constant as long as the orbit is a 
geodesic. E and L represent the energy and the (azimuthal) angular momentum 
per unit mass, respectively. Q is called the Carter constant. Denoting the mass of 
a particle by /i, the energy, the angular momentum and the Carter constant of a 
particle are E = jiE, L = fiL and Q = fJp'Q, respectively. Another notation for the 
Carter constant defined by 

C = Q-{aE- Lf, (2-8) 

is also convenient since C vanishes for orbits in the equatorial plane. We also use 
C = C/n^. 

We can specify an orbit of a particle by using three constants of motion, the 
total energy, angular momentum and Carter constant. Introducing a new parameter 
Xhy dX = dr/U, the equations of motion are given as 

^ = -a(a^ sin^ 9, - L) + !i±^P(r,), (2-9) 

R{r,), (2-10) 

&{cos9^), (2-11) 



dX 
d cos 9z 



dX 



where 



P{r) 
R{r) 
0{cos9) 



^ = -(«^--4^ (2-12) 
dX \ sin^Oz j A ^ ' ^ ' 



Eir"^ + a^)-aL, (2-13) 
[P{r)f - A[r^ + {aE ~ Lf + C], (2-14) 
C - {C + a^{l- E'^) + L^) cos^ e + a^{l- E^) cos^ 9. (2-15) 



It should be noted that the equations for and 9^ are completely decoupled by using 
A. Moreover, R{r) and ©(cos^) are quartic functions of r and cos^, respectively. 

We first consider the radial component of the geodesic equations. When the 
radial motion is bounded by the minimal and the maximal radii Tmin and rmax, ^z(A) 
becomes a periodic function which satisfies rz{X + Ar) = rz{X) with period 

A = 2 -j== (2-16) 

Therefore, we can expand the radial motion in a Fourier series as 

rz{\) = Y.^ne-'''''^^, (2-17) 



where 



fir = 2Tr /Ar. 



(2-18) 



The adiabatic evolution of orbital parameters in the Kerr spacetime 



5 



We can deal with the motion in ^-direction in a similar manner. When the 
minimum of 9 is given by 6'min(< 7r/2), the maximum is ^max = tt— ^min because of the 
symmetry with respect to the equatorial plane. As in the case of the radial motion, 
cos6z{X) becomes a periodic function which satisfies cos 6*2 (A + Ag) = cos (A) with 
period 

d{cos 0) 



cos 6»min 



^0(cos0) 

We can expand cos6z{X) in a Fourier series as 

cos02(A) = ^z„e-*"^«\ 



(2-19) 



(2-20) 



where f2g = 27:/ Ag. 

Next, we consider the t- and 99-components of geodesic equations. Eqs. H2-9() 
and (|2-12() can be integrated as 



tz{\) = t^'\\)+t^'\\) + 
V9.(A) = v5M(A) + <^(^)(A) + 



dtz 
dX 

dip 



A, 



dX 



where 



t^'\X) 
t(^)(A) 



dX 



{rl + a^)P{rz) / (r^ + a^)P{rz) 



A{r, 



Air,, 



dX 



dX 
dX 



c?E sin^ Qz — o,L — (a^E siv? 6^ — aL 

aPjrz) _ I aPjrz 
A{rz) \ A{rz 

L ^ I L 



aE 



sm 



sin^ 6z 



aE 



{■ ■ ■) represents the time average along the geodesic: 



(F(A)) := hm ^ / dX' F{X') 



(2-21) 
(2-22) 



Here, t^^\X) and ip^^\X) are periodic functions with period Ar, while t^^\X) and 
ip^^^X) are those with period Ag. 

§3. The Time Evolution of the Constants of motion 

If the timescale of the orbital evolution due to gravitational radiation reaction is 
much longer than the typical dynamical timescale, we may be able to approximate the 
particle's motion by the geodesic in the background spacetime that is momentarily 
tangential to the orbit (osculating geodesic approximation). Under this assumption, 



6 



N. Sago, T. Tanaka, W. Hikida, K. Ganz, H. Nakano 



we evaluate the change rates of the constants of motion at each moment. For bound 
orbits we can express the change rates of the constants of motion, P = {E, L, Q}, as 



jiinr,ng) [_i^nrOr + n0f2e)A] . (3-1) 



(n,,re)7^(0,0) 



The first term on the right hand side is a time-independent dissipative contribution 
due to radiation reaction, while the others are oscillating. Integrating over a long 
period, the first term becomes dominant. In the same spirit in Ref. 10), here we 
define the 'adiabatic' evolution as an approximation which takes account of only the 
first term. Namely, the adiabatic evolution is solely determined by the time averaged 
change rates of the constants of motion. 

Owing to the argument given in Ref. 10), we can evaluate the averaged change 
rates of the constants of motion by using the radiative field of the metric perturbation 



dP 



lim — - 



dXE—r[hi'"^' 



-T 



(3-2) 



where /ij]^ is the radiative part of the metric perturbation defined by half retarded 



field minus half advanced field, i.e., h^^f^ := {h'^^^ — h^)/2. Radiative field is a 
solution of source-free vacuum Einstein equation. The singular parts contained in 
both retarded and advanced fields cancel out. Therefore we can avoid the tedious 
issue of regularizing the self-force. is a differential operator. 



(3-3) 



This operator with its index lowered reduces to 



•■fiui 



(dah-ys) v?u^ — — {hoi-jU') 



2 ^ " '"^ dT 

ignoring the second order terms. 

3.1. Calculation of dE / dt anddL/dt 
From Eq. H3-2|) . we obtain 



dE 



dXU 



lim 

T^oo 2T 

lim — ^ / dX 



T 



rid 



2T 



lim — ^ 

T^oo 2T 



dX 



E 
'2 



rad 



u<u^ 



d_ 
dX 



'dT 



(/i-^5^xV) +0(^2)^ (3.4) 



h 



E d 
^ '2dX 



(3-5) 



In the last equality, the total derivative terms are neglected. 
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Air) 



sm ( 



(3-6) 



This vector field is a natural extension of the four-velocity of a particle. In fact, it 
satisfies ^^(^(A)) = u^{\). depends only on r and 6. Furthermore, since Uj. and 
uq depend only on r and 9, respectively, we have the relation u^.^^ = Uy^^. 
Using this vector field, we can rewrite Eq. (|3-5|) as 



dE\ -fi 
—— ) = lim — — 

dX / T^oo 2T 



T 



dt ( ^h'^fifu^ 



X'— ►2(A) 



(3-7) 



where we used the fact that E and are independent of t (and ip). 

As shown in Appendix^ (Eq. (|A-54|) ). the radiative field of metric perturbation 
is given by 



1 



2iuj^ 



n: 

^ I ATdown|2 rrdown 



out 1 2 T-rout 



nxU^) I dX UsnTU^{x))u^u^ 



(c.c), (3- 



where A = {Imuj}, k = to — ma/2Mrj^ and r+ = M + \/ M'^ — a^. sU^^^^ix) and 



sn^^°^^\x) are the out-going and down-going mode solutions for /i^j^^, respectively. 
A'^^"* and Nf°^^ are normalization factors, given by Eqs. (|A-52|) and HA-53|) . A bar 
represents complex conjugation. Using this formula, we obtain 



im 



to 



+ (c.c), (3-9) 



where 



^(out/down) := Ari-^'/<i°-'^)r(x),77i°;y''°""^(x)u^(x)n^ 



X . 



For a bound orbit, we can expand 0^ in a Fourier series as: 



,(out/down),^,,,, - -/ — ^Snr^('-)exp 



'(^(A)) 



27r \ dA 



nr,ng 



'./dt, 

'\-dx 



mrirng 



where 



(J. 



mrirng 



dX 



-1 



d^Pz 



m { -j^ ) + Hrfir + neQe ) . 



(3-10) 

)A 
(3-11) 

(3-12) 
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Substituting Eqs. (jHUl) and ^T\l into we obtain: 



1 2 ^mrirng 



^down 



(3-13) 



where fe 



ma 



/2Mr+, {F{t))t ■= liniT-^oo ^ /5rt^*^(0> and 



, (out /down) ^ (out /down)/ -, 

'imnrng = ^ imnrUg {^mrirng)- 



(3-14) 



In a similar manner, the formula for the loss rate of the angular momentum is given 

by 



m 



dL 

dt /. ^ ^ 47ra;3 

* imnrng mnrng 



E 



yout 

^ Imn^ng 



1 2 ^mUrng 



down 



(3-15) 



3.2. Calculation of dQ/dt 

To obtain the change rate of the Carter constant, we need to evaluate 



dQ 



!^ r 

T^OO 2T J_rp 



dX2UK^u^Uhlf]. 



(3-16) 



Using the vector field u'^{x), which was introduced in the previous subsection, we 
obtain 

2K^u^U = lim iK^uf^daih^su^u^) + 2h^su^u'{Kf^.^u'' - K^u%)] , (3-17) 

to the first order in perturbation, excluding total derivative terms with respect to 
r. Those total derivative terms do not contribute after taking a long-time average. 
Furthermore, one can show that the second term also vanishes by using Ki^^p.^-j = 
and Ua;f3 = Ufs-a- After all, we find 



dQ 
dX 



lim 



2 rT 



T^oo 2T J_T 

2 f-T 

lim — — 

T^oo T 

P(r) 



dX 
dX 



2UK^u'^da 



>2(A) 



X 



A 



(3-18) 



To obtain the last term in the last line, the term with m''9^ was rewritten into 
E^^d/dX, and integration by parts was applied. 

Substituting Eqs. 1)3 -Bl) and H3-ll() into Eq. H3-18|) . we obtain: 

■dQ^ —3 rT 



dX 



lim — ^ / dX [ duj 



iniUrng 



UJ — UJ> 



mUrUg j 



imnrng 
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+ 



k 

+(c.c.). 



(3-19) 



xout 



and 



Now we focus on the r-derivative term in the curly brackets. Since 
depend on t and (p only through an exponential function e~*^*+*''"'^, we can write 

^mnrTlg ) 

= f (rziX), COS 6 z{X))6{uJ - tOmnrng) 

' dt 



/.down 



X exp 



mrirng 



+im 



d\ 



X + t^'-\X) + t'^^\X) 



/(r2(A),COS 9z{X))6{uJ - tOmnrne) 



X exp 



irir^rX - iuJmnrngt'-'^\X) + inHf^''^ {X) 



-ineOgX - iuJmnrngt^^\X) + irrup'^^\X) 



(3-20) 



where f{r,cos9) represents the dependence on r and cos0 in (/>yi(x). rz{X), t^^\X) 
and f^^\X) are periodic functions with period Ar, while 9z{X), t^^\X) and (f^^\X) 
are those with period Ag. We introduce two different time variables A,, and Xg. We 
use them instead of A for functions with period A^ and Aq. Then, by using these 
new variables, we can replace the infinitely long time average with a double integral 
over a finite region: 



lim 



1 



T^OO 2T J_rp 

1 



dv 

dX5{uJ - u;mnrng)-T7-dr4'A{ziX)) 



dX 



ArA, 



r^ie Jo 



dr. 



dXr I dXe6{u) - LVmnrne)-^dr\ f{rz{Xr), COS 9 z{Xe)) 



X exp 



iUrQrXr - i^mnrngt'''^'\Xr) + irUif^''' ' (Xr) 



-inefieXe - i(^mnrngt^^\Xe) + imip^'^'iXe) 



(3-21) 



We only need to integrate over one cycle for each of A^ and Xg. Using the relation 



_d_ 

dXr 



^f{rz{Xr),cos9z{X0))eyiY>[-inrOrXr - iuJmnrngt^'"^ +imv3('')]| 



dt^''^ ^ drz ^ d^^'-^ ^ 

xf{rz{Xr),COs6z{Xg))eX.p [-inrf2rXr - iuJmnrngt^''^ + im^j^'')] 

A^-integral in l|3-21|) can be rewritten as 

drz 

dXr-^dr{f{r^{Xr),cos 0z{Xe)) 
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X exp [ - iUrilrK - ii^mrirngt^^' + imif^^>] } 



din I ^ f ^ T 



x{f {rz{Xr), COS 6z{Xr)) exp [ - irirfirK - i^mnrnet'^^^ + irruf'^'^'^]^ }. (3-22) 
Eliminating the r-derivative term from (|3-19() by using the above relations, we obtain 

lim — — - / dX duj y 

imrirn- 



dQ 
dX 



[UJ — LJ- 



mrirng 



yont 

imrirng 



+(c.c.) 



-2/i^ 



A 
A 
1 



dX I ^ 47ra;2 

IranrUg mnrng 



dt + 
dt + 



aP 



A '"^ 
aP 
~A 



d^+inrQr\(t>7 {x) 



x— >2(A) 



dtz_ \ \ ^ 



i70Ut 

(.mrirng 



m 



UJ. 



mrirng 



aP 
~A 



+ 



UJ, 



mnrng 



+ 



UJ. 



mnrng i r^down 



I, i — emnrng\ 

i^mnrng 



(3-23) 



Here we used Eqs. (|3-13p and (|3-15|) in the last equality. Finally, we obtain: 



dQ\ 
dt I 



2^1 



(r2 + a2)P 



A 



dE 



dt \ 



2fj. 



aP 
~A 



3^ - 

imnrn^'^'^^"'r,ng 



z. 



tmnrng I 



+ 



dL 

mnrUg 1 ,7down 



7 \—imnrng\ 
i^mnrUg 



(3-24) 



3.3. Consistency of our formulae in simple cases 

In this subsection, we examine our formulae in a few simple cases. First, we con- 
sider circular orbits. We know that a circular orbit remains circular under radiation 
reaction. This condition fixes dQ/dt for circular orbits as 

dQ _ 2fM{r'^ + a'^)P dE 2fiaP dL 
'dt ~ 



A 



dt 



A dt 



(3-25) 



Since we have Z'^^n^n^^ ~ ^ 7^ in the case of a circular orbit, the last term 

in Eq. H3-24() vanishes. Thus Eq. ()3-24|) is consistent with the above condition that 
a circular orbit remains circular. 

Next, we consider orbits in the equatorial plane. An orbit in the equatorial plane 
should not leave the plane by symmetry. This can be confirmed by rewriting the 
above formula in terms of C. From the definition of uj. 



mnrng 



(|3-12|1 . we obtain the 



following identity: 

2 ^ ^ Tljf 



E 



Attuj 



yout 

(.mnrng 



|2 



mnrng 



^mnrng 



7' 



down 
(mnrng 
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E 



\ d\ 



mrirng 



dtz 
dX 



I I yOUt I 



dE\/dip, 
dt \ dX 



mrirng 



k 



+ 

"(ft 



^ Imnrflg I 



7' 



down 
Imrirng 



2 



+ 



mrirng 



mrirng 



7down 
■'imnrUg 



where we used the the expressions of {dE/dt)t and {dL/dt)t given in Eqs. (|3-13p and 
(|3-15|) . Using this identity, we have 



d£\ _ /dQ 
'dt Xlt 



dt /, 



dt 



3 ne^e 



i,m,nr,ng^^"^"--'^i> 



imnrng \ 



+ 



mUrng 



mnrng 



ydawn 

ttnUrng 



(3-26) 



where we have used the following relations: 

'^') = -a{aE-L) + {a'Ecos'ez) + 

dX 



r1 + 
A 



P 



aE + L+ (Lcoi^ez) + 



aP 



From this equation, it is found that {dC/dt)t = when 9 = 7r/2. Note that we have 



out / down 
imnrUg 



7^ only for = in the case of equatorial orbits. 



§4. Application of our formulation to orbits with small eccentricity and 

inclination 



In this section, as an application of our formulation, we consider a slightly eccen- 
tric orbit with small inclination from the equatorial plane. Since, in this case, we can 
expand an orbit with respect to the eccentricity and inclination, we can analytically 
calculate the change rates of the constants of motion. 

4.1. Orbits 

Here we define ro so that the potential in r-direction i?(r) takes its minimum at 
r = ro: 

dR 
dr 



0. (4-1) 



r=ro 

We denote the outer turning point by ro(l + e). Namely, 



R{ro{l + e)) = 0, 



(4-2) 
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which gives the definition of the eccentricity e. We also define a parameter y = C/L'^, 
which is related to the inclination angle. For orbits in the equatorial plane, we have 
y = 0. Further, we introduce a new parameter v = M/r^. For circular orbits 
V represents the orbital velocity at the Newtonian order. Hence, we regard w as a 
parameter whose power indicates twice the post-Newtonian (PN) order. 

Solving (|4-H) and (|4-2j) for E and L, they are expressed in terms of e and y as 



E = l- -v^ + -V 
2 



2 I ^„,4 „„,5 



qv 



-v' - -y + Iqv" + -qv'^y + qv^e'y, (4-3) 



L = tqv 



63 



2 o 2 

+ (^-l + ^v^-6qv^ + ^v^+'Lq2y^ ^ 

+ ( - ^ - 1^^' + 3qv^ - - \<fv^ + Y'?^^' )y 



-qv e 



13 

+ ( — + Sg?^^ - —v'^ — T-Q.^'v'^ + —qv''']e''y 



n 4 19 2 4 ,63 5N 2 



(4-4) 



where q := a/M. Hereafter we keep terms up to Oiv^e^y) relative to the leading 
order. 

With the initial condition set to rz{X = 0) = ro(l + e), the solution for ^^(A) is 
obtained in an expansion with respect to e as 



=2^(2)1 



(4-5) 



where 



cos f2r\, 



P^'>{1 - cosf^rA) +p(^)(l - cos2r?^A) 



f2j. = vqv 



45 



1 — -w^ + 3qv'^ v'^ q'^v'^ H qv 



33 



1 + -v^ - 6qV^ + 



165 9 
— + 



2 \„,4 



-qv^ 



-qV" - 2q'^v'^ + -^qv'' ]y - ( SqV 



33 



165 

27 2 4 165 
-jQ V + —qv-]e-y 



5 \ „2„ 



(1) 

(2) 



-l-v'^ + 2qv^ - 6v'^ - q^v^ + 20qv^ - {qv^ - 2q^v^ + IQqv^) y, 



We also compute cos 02(A) in a series expansion in y as 

cos0.(A) = Vy[4°)(A) + yc«(A)], 

where 



(4-6) 



.(0) 



.(1) 



1 - -o'l;" 
2" 

1 13 



-q v e ] sin f^gX, 
39 



16 



16 



+ —q'^v^ + —q^v^e^ ] sin QqX + ( —q^v^ + —q^v^e^ ] sin 3^7^ A 



16 



16 
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rov 



1 + -.^. 



3qy3 + ^ 3^2^4 _ ]^^^5 



+ 



6qv'^ + 



81 4 9 4 63 

— V H — q V 

8 2^ 2 

9 2 4 63 5\ 2 

V* + —qv''je'y 



Iq^v^ + ^qv^^y+ (Sqv^ 



Here the solution satisfies the condition, cos02(A = 0) = 0. 

Substituting and cos9z into Eas. H2-9() . H2-12() and ()2-23|) . we obtain 



, dX 



(r) 



. dX 



(0) 



V 



|(2 + 4t;2 - 6gw^ + 17v^ + 3g^f^ - 54qv^) 

+(2 + 6t;2 - Wqv^ + 33v^ + 5g\*^ - W8qv^)e 
+ {3qv^ - + 27qv^)y 

+{5qv^ - 8g^f ^ + 54:qv^)ey^ sin A 

r/3 To 13 s 81 4 13 n 4 135 .\ 



-'^3 ,3 ^„2„,4 



+ ( — 9^^" - 2^"'^" + — 9^^^)62/1 sin 21?^ A 



q^v^r^y 



+ 



1 1 2 
4 2 

1 11 

■ T + -r' 



3 o 5 o 4 r 

-gv H — g w + qv 

4 8 



3g^3 + ^ ^^2^4 ^ ^g^;5^e2} sin2f2,A 



^0 



27 



1 + -w H 1> , , . , 

2 8 " V2 4 "16 

Iq/I 3c\ / „Q „ O A 123 



21 



315 



123 



(4-7) 



(4-8) 



5 \ ^2 



+ (2^^^^ + 2^^^)^ + ( - 3?^^^ + 6g^t^^ - —qv^y^y 

I (-2 + 2gu - lOv^ + 18g7;^) + (-2 + 2qv - 12v^ + 2Aqv^)e 

— {qv + 9qv^)y — {qv + 12qv^)ey^ sin i?^ A 
+{( - ^9^^ + ^^'^ - ^^'"^) ^ + (^9^ + ^9^^) sin 2f2,.A 



(4-9) 



y 

rov 



^qv + -V 

1,32 4\,9242 

- + -qv)+-qve 



sin 2f2gX, 



(4-10) 
(4-11) 



2 „„,3 I 2''', ,4 9 5 



1 H — v^ — qv^ H V 

2 8 



2^ 



81 



27 



5 \„2 



1 v^ + 2qv'' v^ H qv° e 



15 



(4-12) 



4.2. Calculation of 



In order to obtain the averaged change rates of the energy, angular momentum 

•out /dow: 
(.mrirng 



and Carter constant, we have to calculate Z'^^l'^^'^ defined by Eq. (|3-14|) with 
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Eq. (|3-11|) . Integrating Eq. (|3-10|) with respect to A, we obtain 



A(out/down) 



= / dX(j) 



(out/down) 



(^(A)) 



^(out/down) / ^4^yZ^^^{out/down)^^^ / 



,T (out /down) 



where 



T"^(x) =fl [ d-r ^ 7/"7/^W(^ _ z{t)). 

J \/-g{x) 



(4-13) 
(4-14) 



is the energy momentum tensor of a mono-pole particle of mass /i. Using the relation 
given m can be also expressed in the familiar form which 

appears as an integration over the source term in the standard Teukolsky formahsm 
as 



,^(out/down) 
^ (out /down) -''s 



Cs I d^x./^^ sRf^''''\r)sZA{9,'p)e''^'sf{x), (4-15) 



where sT{x) is a projected energy momentum tensor defined by gT := sT^^T^^ 

with (|A-9p . and sR^A^^^'^\'r){= -sR^^^^^'^°'"'^\r)) and are, respectively, the 

radial mode functions and the spheroidal harmonics introduced in Appendix 1X1 2. 

In the following discussion we concentrate on the case with s = —2. Substituting 
the explicit forms of the energy momentum tensor and the projection operator -2T^u, 
we obtain 



/oo 
dte 
-oo 



iu^t-imv{t)j(m/np) ^^^^^ ^ ^ ^ ^ 



with 



-^rhnO 



+ A 



mmO ) 



27rZ\2 



dRA . . A. . \ d'^RA . 

CnnZ^zC\[z^4{z-^SA)}, 



r{t),e=e{t) 



r- 7-^ 



Z Z Crnm Sa 



^ + + z-^^cISa - ^{z-^ - z-^)asmeSA 



2tt 



.(K\ 2iK 



Z^Cn 



£2'S'yi + iasin^(2: — z 
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2tt 
1 



Cn 



2tt 



Z Z Crt 



iSa (^i^ + z ^ 



where Sa represents -2Sa{Q) defined in Appendix. IXl and 



z 
K 



r + ia cos o 

de + - 
sm a 



- ma, 

auj sin 6 + s cot ( 



Here dagger (f) means an operation that transforms (m, uj) to (— m, — w). The radial 
functions and the spheroidal harmonics appearing in the above equations can be 
evaluated analytically, as shown in Appendices El ^-nd [O For a bound orbit, since 



e-i^v{t)xi^^/''P)(^r{t),e{t)) is a double periodic function, Z 



^(out/down) 



spectrum as 



^out/down 



\ ^ „out/down r-f 
nr,ne 



U). 



mrirng ) ) 



has a discrete 
(4-17) 



where the coefficients Z\ 



•out / down 
imnrllg 



are those already introduced in Eq. (|l-{-14j) with 
Eq. Although we cannot show all the processes explicitly here, it is straight 

forward to calculate Z'^^^^^™ for each uJmnrne by substituting the analytic expan- 
sions of the orbits, the radial functions and the spheroidal harmonics. 

4.3. Results 

Substituting Z'^^^'^°'^^ obtained by following the scheme explained in the pre- 
ceding subsection into Eqs. ()3-13() . ()3-15() and 1)3 ■24|) . we obtain: 



'dt 



32 

y 

1 



MJ 
1247 



10 



336 
/ 44711 

~ V 9072 

r 277 4001 

n 



73 
12' 



47rJ 
/3749 



33 2 \ 4 I 
16 / V 336 



8191 



672 



24 
+ (^2q 
73 



5 /3583 

84 V 48 

2 1091291 \ 4 



9072 / 



457 

/ 58487 
V 672 



364337 
1344 



+ 1^9^ - 
/457 3 



527 



2„,4 



3749 



-q V 



96 

5407 , 4 
48 ^ 



-qv y 



672 

58487 2 

qv e V 

1344 ^ ' ^ 



(4-18) 



dL 
'dt 



32 



M 
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1247 



—q-4Tr]v' 



+ 



33 2 \ 4 I 

Q ]V + 

9072 16^ ' 



336 
^44711 

r 51 _ 17203 

\Y ~ 672 

/929 2 1680185 \ . 
V 32 ^ 18144 y 



/ 417 _ 8191 
V'56"^~ "672 
781 369 

T2-^ + ^")^ 



/1809 



224 



■q - 



48373 
336 



TT ]V 



1 1247 2 
+ ^-2 + W^ + 



61 



-g - 27r t;' 



213 
"32"' 



44711 \ 



18144 



51 17203 2 /1513 



/4301 
V 224 ' 



8191 

1344' 



TT \V 



'^^ 16 1344 



369 



+ 



/ 1680185 5981 



V 16 ^-T^"'" 



dQ\ 64 3 , 



dt/t 



36288 



743 



64 



( 168g - 



48373 
672 



TT \v >e y 



(4-19) 



1 ■ 
+ 

+< -5 ri'" 



/1637 



- qv V 

^ 336 

/439 2 _ 129193 

Vls"^ " 18144 ~ 

43 51 2425 



-g — 47r u 



336 

- 47rg ) ti'* + 
2 



/ 151765 4159 



+ 



/ 453601 

V 4536 
/ 141049 



224 
3631 



V 18144 

/ 14869 337 

q — TT 1 V 



224 



672 
3 



-TT 



33 
16' 



32 

38029 



9072 " 

1 ^ 1637 , 
+ < -qv H ;;—qv 



672 



-TT 



369 

—Trq ]v 
929 3 \ 5 



2' 672 
_/ 151765 



/1355 



96 



— 27rg I v 



V 36288 ^ 32 ^ 



+ 



[51 

ll6^"+ 448 
/ 141049 



14869 3 /369 

■qv +(^vrg. 



33257 2 



5981 



18144 " 64 " 
Prom the above results we can compute 



V" >e y 



(4-20) 



dQ\ / 2fj,{r- + a^)P \ /dE\ / 2fiaP \ I dL 



t \ 
64 3 fi 2 
5 



+ 



37 13435 2 

\ - 

6 672 

/625117 337 



/1561 



335 



-TT 



12096 



32 



f 46827 1355 

■TT — 



V 448 



672 



q ]v" 
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335 
l6 



■qv 



7559 
l92 



,2„.4 



■q V 



1355 
1344 



qv \y 



(4-21) 



The left hand side of this equation vanishes for circular orbits. ^^-^ In fact, the right 
hand side vanishes for e = 0. When we did not know how to compute {dQ/dt)t, 
the best guess that we could do for {dQ/dt)t was to assume that the left hand 
side vanishes for general orbits. ^^-^ Therefore this combination represents the errors 
coming from this hand- waving working hypothesis. Wc can also compute 



{dt\ 



dQ 
dt 



- 2{aE - L) 



64 3 6 
5 



743 2 

v'^ - 

336 

/ 129193 
I 18144 



dE\ _ 

85 

8 ^ 
307 

"96 



dL 
It 



■q 



US 2425 2 /337 

-l--^ V + TT - 

\ 8 224 \ 8 
/ 453601 7849 



/2553 4159 
1793 



16 



q ]v 



672 

3 



TT \V 



+ 



V 4536 
/3421 



224 



q 



192 
38029 
672 



V 



TT ]V 



(4-22) 



Since y = (i.e., C = 0) corresponds to = 7r/2, C does not evolve for equatorial 
orbits. This is consistent with the requirement from the symmetry that an orbit in 

the equatorial plane stays in the equatorial plane. 

Finally, we consider the evolution of an inclination angle l defined by,-*^^^ 



L 



cos L 



(4-23) 



which, roughly speaking, represents an angle between the normal vector of an orbital 
plane and the rotational axis of the central black hole, but this is not the unique 
definition. Although the definition of inclination angle can be changed at will to 
some extent in Kerr case, thus defined inclination angle reduces correctly to the 
usual one in the q = Schwarzschild limit. Taking the average of the time derivative 
of cos L, we obtain 



dcos L 
dt 



2(L2 + C)i 

32fi^v^ 
5L2(l+y)l 



dL 
'dt 



C-L 



qy 



dC_\ 
13 . 1779 



61 3 4 
24 96^ 



431 

Te"' 



775 
192 



qv 



224 
22431 
' 224 



(4-24) 



Substituting g = into this equation, we can confirm that l does not change in the 
case of Schwarzschild limit, which must be so because of the spherical symmetry of 
Schwarzschild spacetime. 
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§5. Summary 

In this paper, we have considered a scheme to evaluate the change rates of 
the orbital parameters of a particle orbiting Kerr black hole under the adiabatic 
approximation. We have adopted the method proposed by Mino/'^) in which we use 
the radiative field instead of the retarded field in order to compute the change rates 
for "the constants of motion" due to radiation reaction approximately. Based on 
Mino's method, we have developed a simplified scheme to evaluate the long term 
average of the change rates. Applying our new scheme, we have performed explicit 
calculations to present analytic formulas of change rates, {dE/dt)t, {dL/dt)t and 
{dQ/dt)t, for orbits with small eccentricity and inclination angle. 

Here we used the expansions with respect to the post-Newtonian order, the 
eccentricity and the inclination angle in evaluating {dE/dt)t, {dL/dt)t and {dQ/dt)t- 
As a next step therefore we need to examine how large parameter region is covered 
by our formulae with a sufficient accuracy. As for the inclination, we recently found 
a formulation to obtain the analytic formulae for the change rates without assuming 
a small inclination angle. On the other hand, it is almost certain that we need 
numerical calculation for the cases with a large eccentricity. Drasco and Hughes^^^ 
developed a numerical code to calculate the gravitational wave fluxes of energy and 
azimuthal angular momentum evaluated at infinity and at the event horizon for 
general geodesic orbits. Fujita and Tagoshi also developed a numerical code based 
on an analytic method of solving the radial Teukolsky equation. By applying such 
codes to our scheme, we can evaluate the time-averaged change rate of the Carter 
constant for general orbits, although computational cost will not be small because 
we need to take into account a large number of frequency modes. 

Once we obtain the change rates of "the constants of motion", as a next step, 
we want to use them to trace the evolution of orbits. Some strategies to solve the 
orbital evolution taking into account the radiation reaction effects were proposed in 
Refs. 17) and 18). However, it should be noted that the adiabatic approximation 
used here contains only the dissipative part of the self-force on a particle, and it does 
not contain the conservative part. In general, the conservative part also contributes 
to the secular evolution of orbits, thought it is not the dominant part in the limit 
/X — 0. Therefore the adiabatic approximation may not be sufficient to evaluate the 
orbital evolution. 

Recently, Pound, Poisson and Nickel showed that the conservative part of the 
self-force can produce significant shifts in orbital phases in an analogous problem 
with a charged particle in electromagnetism.^^^ They suggested that the conser- 
vative contribution to the phase shift is relatively large in weak field, slow motion 
cases, while it is suppressed in strong field, rapid motion cases. Furthermore, there 
are different types of effects higher order in /x which may produce significant shifts 
in phases. Therefore it is important to quantify the range of validity of the adiabatic 
approximation for appropriate applications of the results obtained in this paper. 
Although it requires computing second order perturbations in fi in order to under- 
stand the whole effects which potentially give phase shifts greater than 0(1), some 
of effects can be evaluated by studying the first order self-force at each moment 
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without averaging over a long period. We will come back to this issue in one of our 
forthcoming papers. 
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Appendix A 

Radiative solution for the metric perturbation 



In this Appendix, we give a brief review on Teukolsky formalism, followed by a 
derivation of the radiative Green function of the linearized Einstein equations. This 
derivation is based on Refs. 22), 23) and 8). 

A.l. Teukolsky equation 

As a master variable we consider the Teukolsky functions defined by 

-Cai3'ysl''rn'^l^m^, s = 2, 
—z^Ca/3'yS'n"m^n'^7fi^, s = —2, 



(A-1) 



where 



1 

2z 



1 

2z 



1 
2 



z 



4^/2\ Az^ z^ ^ A 



,(A-2) 



z •— r + ia cos 0, A := r'^ — 2Mr + a^, and U := r'^ + c? cos^ Q. P„ and Lg are the 
differential operators defined by 



Vr, := dr ^ 7 — -Ot + + 



A 



2n[r - M) 



do 



sm( 



— ia sin Qdt -|- s cot ^, 



(A-3) 
(A-4) 



and a dagger (f) acting on an operator means transformation of (Of, ^^p) {—dt, —5^), 
which reduces to the one defined in the main text by {uj,m) — (— w, — m) under the 
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assumption of Fourier expansion. The Teukolsky functions satisfy a separable partial 
differential equation'^^-* 

,0 = ^ttU sf, 



where 



and sO is the Teukolsky differential operator, 

sO := sOr + sOe, 

with 

^2 4Mar 



s^r 



-2s ( '- -r]dt + s 



d-t + A-^dr{A^+^dr)--dl 



A 



dtd^ + 



(A-5) 
(A-6) 

(A-7) 



2sa{r - M) 
A ^ 



.Oe:=a' sin^ 9^ + J- de (sin ede) + -^dl + 

sm u sm^ 9 ^ sm^ 9 

—2isa cos 9 dt — cot^ 0, 



and 



V2 



1 z 



(A-9) 



We consider the following forms of expansions for and ^T: 



oo 



/oo 



27r 



where A := {/ma;}. Substituting these expressions into the Teukolsky equation 
HA-5|1 . we obtain equations separated for the radial and angular parts as 



A Z\'*+^— H h ^isLor - A 

dr \ dr I A 



,XA{r) = sr^,(A-10) 



1 d f . ^ d\ 22-2/1 {m + s cos 
sm6'— - a^uj'^ sm^ 9 - ^ 



sin^d^V d9 



2 a\2 



sm 
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-2aujs cos + X + s + 2amuj 



,5^(0) =0, (A-11) 



where K := (r^ + a^)u; — ma and A := sEem{o,^) ~ s{s + 1) + o^u;^ — 2amuj. The 
eigenvalue sE^micL^) is determined by solving Ea. HA-llI) as an eigenvalue problem 
imposing regular boundary conditions on sSy[{9) at 9 = ib7r/2. Here i is an index 
that labels difFernt eigen values. We also give a brief review on how to solve this 
equation analytically in Appendix C. 

A. 2. Mode functions 

We write mode functions for the Teukolsky equation HA-5|) in the form 

s^A := sRA{r) sZA{e, (^)e-^-*, (A-12) 

where <ji?yi(r) is a homogeneous solution of the radial Teukolsky equation ()A-10|) . 
and sZA{0,ip) is the spheroidal harmonics 

sZAie, V) = ^ sSA{e)e'"''', (A-13) 
V 27r 

normalized as 

/ desmeisSAiO)]"^ = 1. (A-14) 

Using the symmetry of the radial equation (|A-10() under the simultaneous operations 
of the complex conjugation and the transformation of (m, u) — > {—m, —uj), we impose 

sRa= sR\, (A-15) 

where a dagger (f) acting on a mode function means transformation of {uj,m) — > 
(—uj,—m). In a similar manner, by virtue of the symmetries of Eq. HA-ll|) . we 
arrange the spheroidal harmonics to satisfy 

sZa = {-ir -sZ\. (A-16) 

In our later discussions, we also need the well-known Teukolsky-Starobinsky identi- 
ties: 

.sRa = sUsRa, (for \s\ = 2), (A-17) 

with 

~2U := ^Vl ,U := ^A'T>1'^\ (A-IS) 

where 

C = [((A + s{s + l)f + 4acjm - ^a^oj'^){{X + s(s + 1)- 2f + maum - SGa^w^} 
+{2\ + 2s{s + 1)- l)(96a^a;2 - AMujm) - lUa^uj^]^/^ + 12iu;M, (A-19) 

and A is a factor which depends on how we normalize the radial functions. In this 
paper, we simply adopt A = 1. (This convention is the one used in Ref. 24)). 

Now we discuss how to construct mode functions for metric perturbations from 
mode functions of the Teukolsky equation. The basic idea owes to Chrzanowski.^^) 
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Here we follow a more rigorous approach taken by Wald.^^^ Using the relation (|A-1|) . 
the Teukolsky equation ()A-5p is rewritten as 



(A-20) 



On the other hand, operating ^Tq,^ on the linearlized Einstein equation, which we 
schematically denote as G'^^^^'^h^^ = AnTaf^, we obtain 



T. 



(A-21) 



From the comparison of these equations, we find an identity at the operator level: 

1 



Here we define O*^" , the adjoint of an operator O^^^ , so as to satisfy 



gY^yO*i"'X(rx, 



(A-22) 



(A-23) 



for arbitrary scalar field X and tensor field Y^jj. The definition of the adjoint op- 
erators for different types of tensor operators is a straight forward generalization of 
this definition. It will be worth noting y/—g d^x = sin 6 U dt dr d9 dip, and 



(AB)* = B*A*, VI = -U-^vl^U, C* 
By taking adjoint of each side in Eq. ()A-22|) . we obtain 



(A-24) 



(A-25) 



Here we used the fact that the linearlized Einstein operator G'^^^'^ is self-adjoint, 
i.e., G*'^l^^''' = G'^^^'". Then, from the definition of ^O,. and sOe given in Eqs. (fX^ . 
it is easy to see that 



Therefore we have {IJ gO)* -sRa s^a e"*"^* = 0, whihc means that 



■iuit 



0. 



(A-26) 
(A-27) 



Here the explicit form of the adjoint operators ^r*^ are 



fiu 



zz-^ z 



(A-28) 
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1 z 



^. (A-29) 



Hence, 

sIIa^iiu ■= Cs sT*^s^A, (A-30) 

with 

sf?A = -sRAsZAe-'^', (A-31) 

is a complex- valued homogeneous solution of the linearized Einstein equations. Here 
(g is a numerical coefficient which we determine so as to satisfy 

sD^"' ^{Aa snA,^^^ + Aa sHa,^^) = Yj-^^ 

A A 

for any complex- valued amplitude of each mode, A a- Here the coplex conjugate term 
in parentheses is necessary to make the metric perturbation real. Using Eqs. ()A-2|) 
and dSSni), we can verify 

-2l?'^"-2t;, = 0, al^'^V;. = 2?o', (A-33) 
-2I?^"2t;, = ^£_i£oA/:2, 2l?'^V;. = 0. (A-34) 

In literature sfpi is used to represent what we denote here by st*^. The difference 
arises because we use the notation for the differential operators without assuming 
that they always act on a single Fourier mode. Namely, instead of writing {—iuj, im), 
we are using here {dt,d^). The complex conjugation of the former gives rise a flip 
of signature, while that of the latter does not. With the aid of the above relations 
HA-17|) and HA-34|I . we find that the complex conjugate terms vanish to obtain 

2D^'^2nA,^u = C2 C 2i^yl, 
_2l)^'^-2iI/l,M. = ^^2nA. (A-35) 

lb 

Thus the normalization constants are fixed as 

C2 = i C-2 = f. (A-36) 

A. 3. Radiative field 

Here we explain a method of constructing radiative field for metric perturba- 
tions. Radiative field is a homogeneous solution of field equations. Hence, once we 
obtain the radiative field for the Teukolsky function, it can be easily transformed 
into that for metric perturbations by using the relations established in the preceding 
subsection. We therefore first derive the radiative field for the Teukolsky function. 

The retarded Green function of the Teukolsky function is defined as a solution 

of 

,OsG{x,x') = ^-^-^^^^, (A.37) 
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with the retarded boundary condition: sG{x, x') = for t < t'. We write the retarded 
Green function of the Teukolsky equation in the form of Fourier-harmonic expansion 

as 



gCr(x, X 



duj 
2^ 



sgA{r, r')sZA{0, ip)sZa{9', ip')e~ 



iuj(t-t') 



(A-38) 



Then the radial part of the Green function sgAif^r') is given by 

1 



s9A{r, r') 



Rj{T)sRX{r')e{r-r') 
+sI^2{r)sR7{r')e{r' -r) 



with the Wronskian defined by 



^^A{r)^^sR7{r)-sR7{r)^^sI^A{r) 



(A-39) 



(A-40) 



The advanced Green function can be constructed in a similar manner just by 
replacing "in" and "up" with "out" and "down" , respectively. Then it is easy to show 
that the radiative Green function has a simple structure which does not contain any 
step function 9{r — r'). To show this, let us start with the following expression for 
the radial part of the radiative Green function for r > r': 



,Rj{r),RXir') ,Rf^^{r)sRT{r') 



WisKX^sRj) WisRT^sR'y 



(A.41) 



We rewrite this expression in terms of the down-field and the out-held, eliminating 
sRfir) and ,R°/^{r'){= Z\-"_^i?|^(r')) in Eq. (fFTHl . Hence we expand sRj and 
sRT as 



Tjup _ pout , n pdown 

sRT = 7 sRj + 5 sR^X ■ 



(A-42) 



Taking the Wronskians of both sides of Eqs. ()A-42|) with appropriate radial functions, 
one can easily obtain 



up Tjdown N 



W{sRJ,sRa 



out ndown \ 



aW{sRT,sR' 



W{sRJ,sRa 



30Ut\ 



WisRT^sR'X) =iW{sRj,sR'X), W{sRT,sR7) = SW{sR}X,sR7' 



A ^ 



Substituting these relations, the expression (|A-41() reduces to 

rad / /\ ^ ^(t') 

sQa {r,r') 



2WisR'X, sRjWisRT^sRT''' 



WisRT\..R''' 



A ^s^iA) sRA^'^i'^) ~sRA^'^{r) 

+w{sR7,sRa 



up Ddown^^^out(^^_^^out^^/^ 



(A-43) 
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We can do an analogous reduction for r < r', and the result turns out to be the same 
as that for r > r' . Namely, the step functions which was present in the retarded 
and the advanced Green functions do not appear in the radiative Green function. 
This is consistent with the fact that the radiative field is a source-free homogeneous 
solution. 

Since the radiative field is a homogeneous solution, we can use the method for 
reconstruction of metric perturbation explained in the preceding subsection. When 
we consider the metric perturbation by a point mass, the energy-momentum tensor 
is given by ()4-14|) . In this case it is easy to verify that the radiative field of the metric 
perturbations is given by 



with 



im 

+^down^^down(^) / dr ,i7;}°-^-(z(A))n"n^ |> + (cc), (A.44) 



T-down 



down 



CsWisRJ^, sRf)W{sRT,sRf 



?down \ 



In fact, if we apply sT^^^ ^ we correctly recover s^^^{x) 
A^{x')sT{x')d'^x. To show this, we also used 



(A-45) 



j G''"^{x,x')S{x') 



-QsfiA sTd'^x = CsJ V^[sr;, s^AjT^^'d'^x 

= ^i [ dTu^'u^nA,^,u{z{T)). 



(A-46) 



It is more convenient to rewrite Ng written in terms of Wronskians by using the 
coefficients in the asymptotic forms of radial functions. The radial functions take 
the asymptotic forms. 



riin 

s^A 



sR7 



^trans^-Sg-ifcr* 



+ sBfr-^-'-'e'^'\ 



Dine -1 - 
s^A 



/^trans — 2s— 1 iojr* 



for r* 

for r* 

for r* 

for r* 



oo, 
— oo, 

• oo, 

• — oo. 



(A-47) 
(A-48) 



where r* is the tortoise coordinate defined by dr* /dr = (r^ + a^)/A. Using the rela- 
tions sRT = ^~^-sRT and sRf^'^ = A-'_sR7, we can describe the asym ptotic 
forms of out- and down- fields with the same coefficients that appear in Eqs. HA-47|) 
and ()A-48|) . Then, the Wronskians that we need to evaluate are 



w{sR'X,sRT) = ^^'^sBTsOa 



trans 



w{sRT,sRTn 



Oo, , /^trans oinc 
-ZILO -s'-^A -s^A 1 
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W{sRT,sB}X) = -4ikMr^t,,sBT''-sBT', 
WisRf, sRTn = -2iio ^sCT'"' sCT''', (A-49) 

where Kg '■= 1 — is(r_|_ — M)/2kMr^. The coefficients with (— s)-spin can be erased 
by using the Teukolsky-Starobinsky identities (|A-17|) . Substituting the asymptotic 
forms (|A-47|) and (|A-48|) into Eqs. HA-17|) . we obtain 

C ■ / 1 \^ r 

Dine rjtrans i I rjtrans 
, -on A = —— r 



/-ttrans oref 



Using the above relations, the coefficients Ms are rewritten as 

1 1 

A rout I A rout 1 2 * rdown i Ardown 1 2 A 1 ^ 

with 

23s-2j^2s+2 

1^"'!' = |C|./2-l |^5inc|2' (^-52) 

9-35-2 2,-25+21/^15/2+1 I Dtrans|2 

I A7down|2 ^ ^ ^ I*- 1 U-°yl I ,'A 

' |K2|^/2-i|Ki|^(2Mr+)2^-i Is^rPI^'^r^'P' 

Hence, we finally obtain 

h^;^ = f, j dcoY,^{Nr\nZuix) [ ^4>r{r) 

''em 

where 

^(out/down) := Ari°-t/d°-'^)i7(z(r)),i7j;5/''°"'°)(z(T))n^(r)n^(r), (A-55) 
whose extension to a field is 1;^^"*'^'^°™'^^ (j;) defined in Eq. H3-1U() . 

Appendix B 

Mano-Suzuki- Takasugi method 

Mano, Suzuki and Takasugi formulated a method of constructing a homoge- 
neous solution for the radial Teukolsky equation in two kinds of series by using the 
Coulomb wave function and the hyper geometric functions. ^^-'"^^^ By applying this 
method under slow motion approximation, we can express homogeneous solutions in 
an analytic form. Furthermore, this method determines the asymptotic amplitudes 
of homogeneous solutions without numerical integration. This allows us to compute 
the gravitational wave flux at infinity and on the horizon with a high accuracy. 
We summarize this method in this appendix. 
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B.l. Outer solution of radial Teukolsky equation 

According to 24)-26), we can expand s-R^"' ^ homogeneous solution of the radial 
Teukolsky equation ()A-10|) . in terms of the Coulomb wave functions as 



_ r{,, + l-s + ie) / eKys-ie^ 

XiFi(n + 1/ + 1 - s + ie, 2n + 2z/ + 2; 2iz), (B-1) 

where e = 2Moj, e+ = e + r, r = ^-^(e - ma/M), k = ^/l - {a/Mf, := 
r{x + n)/r{x), z := uj{r — r_), and r_ = M — \/ M'^ — a^. The coefficients 
satisfies the following three term recurrence relation, 

<d + «^ + 7X'!i = 0, (B-2) 

where 

ieK{n + V + 1 + s + ie){n + V + 1 + s — ie){n + V + 1 + ir) 



(n + z/+ l)(2n + 2z/ + 3) 



= -A - s{s + l) + {n + u){n + iy + l)+e^ + e{e- mq) + "^"^^^'^ + ''^ 



(n + i'){n + J/ + 1) ' 

^ ieK(n + — s + ie)(n + 1/ — s — ie)(n + z/ — ir) tr, n\ 

" (n + z.)(2n + 2z.- 1) ' ^ ' 

and g = a/M. The renormalized angular momentum i/ is determined by the condi- 
tions 

an ten a„' ten 

hm n-^ = hm n^;;^ = — — . (B-4) 

n— 1 n+1 



Under this condition, the series of Coulomb wave functions ()B-1|) converges for any 
r > r+. 



From the equations in ()B-3|) . we can show that = 7^ and j3_^ = /?^. 

By using these relations, we can find that "^'^ = a^^ and 

— U—1,S ■ — I/ — l,s 

an ten an ten 
lim n ^ — = , lim n --, — = . B-5 

This fact shows that gRc^ ^^ is also a solution of the radial Teukolsky equation, 
which converges within the region r > r+. 

B.2. In-going and up-going solutions 

The in-going solution of the radial Teukolsky equation is given in terms of the 



Coulomb type solutions (|B-1|) as 

sRX = Ase'^'^iKs,, sRc + Ks,-u-i sRc"~'), (B-6) 
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where 



yen/ -T (— 1 — 2?,e^ 



r{u- l + ie) 



r(i/ + 3 + ie) 



, ^-2 = 1, 



(B-7) 



{2eKy 



(z^ + 1 + iT)r{y + 1 + s + ie)r 

r{l-s- 2ie+)r(n + 2v + 2)r(n + 2v + l) 

r{r + 1/+! - s + ie)r{u + 1- s- ie)r{u + 1 - ir) 
„ (r + 2z/ + + 1 + s + ie)n{iy + 1 + ir 



E(-i) 

,n=r 
r 

E 



(n - r)! (u + 1 - s - ie)„(z> + 1 - ir). 



-1)" 



(r - n)! (r + 2z/ + 2)n (j^ + 1 - s + ie 



(B-8) 



Here r is an arbitrary integer and Kg^i, is independent of the choice of r. 



Next, we consider the up-going solution. gRr: can be divided into two parts as 



(B-9) 



where 



g-7regi7r(i/H-l-s)g-ii 



^ ' V zJ r(u + l + s-i€) 



r{v + 1 + s - ie) 

oo 

^ (2ii)"a;;'^!Z'(n + z^ + 1 - s + ie, 2n + 2zy + 2; 2ii), (B-10) 



X f;(2zl)" ^" + ^ + ^-^^^- a--- 



(z/ + 1 - s + ie) 
xip-(n + z^ + 1 + s - ie, 2n + 2z/ + 2; -2ii), 



(B-ll) 



and •P'(a, c; x) is the irregular confluent hypergeometric function. Prom the asymp- 
totic form of •Z'(a, c; x), 



^{a,c;x)^x — ^ oo), 
the asymptotic forms of gR'^ and gR- become 



where 



41/ -1 -i(2;+eln2;) 



^R'L 



-7re/2 i7r(:/+l-s)/22S-l-i£ -^(^ + 1 - S + ie) ^ 



r(i/ + 1 + s - ie) 

„ (z/ + 1 + s - ie)„ 



(z^ + 1 — s + ie)„ 



(B-12) 
(B-13) 

(B-14) 
(B-15) 
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This shows that sR- {sR+) satisfies the up-going (down-coming) boundary condition 
at infinity. So we can take the up-going solution as 



BssR-, 



where B2 = Clo'^^ and B^2 = 1- Taking the limit r* 
by means of the asymptotic form of r* , 



(B.16) 

±00 in Eqs. (fFfij) and 



u>r 
kr* 



z + e In z — e In e (r 



00 



e+ In(-x) + Ke+ + 



2Ke+ 

1 + K 



In K (r — > r+), 



(B-17) 
(B-18) 



we find that the coefficients which appear in the asymptotic forms of Eqs. HA-471) 
and ()A-48|) are given by 



Dine 



Dtrans 
s^A 



/^trans 



Ase'' 



ijj 



A, 



2s 



le 



, sin7r(z/ — s + ze) 
sin7r(z/ + s — ie) 



n=— 00 
-1— 2s Ae In e 



—00 

A".. 



(B-ID) 

(B-20) 
(B-21) 



Appendix C 

Spheroidal harmonics 



Here, we review the formalism to represent the spin-weighted spheroidal harmon- 
ics in a series of Jacobi polynomials based on Ref. 27), which was slightly improved 
in Ref. 16). 

We first transform the angular part of the Teukolsky equation ()A-11|) as 



(l-x^) 



d 
dx 



vr? + + 2msx 
1 — 



2s^X + sEtmii) 



(C-1) 



where ^ = au^x = cosO and sEimiO = A + s{s + 1) — + 2m^. The angular 
function s5'|^(x) is called the spin-weighted spheroidal harmonics. Equation HC-1|1 
is a Sturm-Liouville type eigenvalue equation with regular boundary conditions at 
X = ±1. Since there are a countable number of eigenvalues for fixed parameters s, m 
and ^, we introduced an index i starting with max(|r?T-|, |s|) as such a label that sorts 
the eigenvalues sEgmiO ™ ascending order. When ^ = 0, sSf^{x) is reduced to 
the spin- weighted spherical harmonics, and the eigenvalue sEimH) becomes t{l + 1). 
We normalize the amplitude of s»S'|^(x) as 



sinOde = 1. 



(C-2) 
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The differential equation ()(M|) has singularities at x = ±1 and at x = oo. We 
transform the angular function as 



and 



1 - x\ 2 /I + x\ 2 



1 - X\ 2 /I + X 



■iV£m{x) , 



(C-3) 



(C-4) 



where a = |m + s| and (3 = \m — s\. Then, Eq. ()(M|) becomes 

(1 - x2) MLix) + [/?-a-(2 + a + /3)x] MU^) 



+ 



a + (3 f a + (3 



+ 1 



sUim(x) 



and 



e [-2(1 - x^) sU'ernix) + (a + /? + 2s + 2)x sUim{x) 
-iC + P - a) sUemix)] , 

(1 - x2) ,y/;^(x) + [p-a-{2 + a + /3)x] ,y/^(x) 



(0-5) 



+ 



a + P ( a + (3 



+ 1 



iiV£m(x) 



= i [2(1 - x^) ,y/^(x) - (a + /3 - 2s + 2)x,y,„(x) 

-(e-/3 + a),y,^(x)] . (C-6) 

From Eqs. (|C-3|) and HC-4|) . we find 

sVimix) = exp(2^x) sUim{x) ■ (C-7) 

When ^ = 0, the right-hand sides of Eqs. (|C-5p and (jC-Sp are zero, and they 
reduce to the differential equation satisfied by the Jacobi polynomials, 

(1 - x2) P^"'^)"(x) + [(3-a-{2 + a + (3)x] Pi"'^)'(x) 

+n(n + a + /3 + 1) P^"'^) (x) = 0. (C-8) 

In this limit, the eigenvalue sEnmiCj in the equation ()C-5|) becomes f(£ + 1), where 
n = £ — [a + (3) /2 = i — max(| m |, | s |). Here, the Jacobi polynomials are defined 
by the Rodrigue's formula by 



-1)" 



2"n! 



(l-x)-"(l+x)- 



d_ 

dx 



(l-x)"+"(l + x)^+"l . (C-9) 



Now, we expand gUemix) and sV£m(a^) in a series of Jacobi polynomials: 

oo 

sUUx) = E ^^t^(0 , (c-10) 

n=0 

oo 

.V^.™(x)= j;.<)p^'3)(x). (C-11) 



n=0 
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The expansion coefficients sA^^{^) and sB^^{^) satisfy the recurrence relations 

«^°^42(e) + /?^°^42(e) = o, (c.12) 

+ sA^lliO + 7^"^ sAf-'\i) = 0, (n > 1) , (C-13) 



with 



4^(n + a + l)(n + /3 + l)(n + (a + /?)/2 + l-s) 



{2n + a + l3 + 2){2n + a + l3 + 2>) 

,3'"':= ,E,„(« + «^-(„ + 2i^) + + i 

2is{a - 0){a + 0) 



+ 



(2n + a + /3)(2n + a + /? + 2) 



and 



„(n)._ ^in{n + a + P){n + {a + 0)/2 + s) 

^ ■ (2n + a + /?- l)(2n + Q + /3) ' ^ ' 



+ + 7^") sBt^'\i) = 0, (n > 1) , (C-15) 



with 



4e(n + a + l)(n + ^ + l)(n + (q + /3)/2 + 1 + s) 



(2n + Q + /3 + 2)(2n + a + /? + 3) 
2^s(a-/3)(a + /?) 



+ 



{2n + a + 0){2n + a + (3 + 2)'' 



-in) ^^n{n + a + 0){n + {a + 0)/2-s) 

^ ■ {2n + a + -l){2n + a + P) ' ^ ' 

The eigenvalues sEimiO ^-^'^ determined in a way similar to the renormalized 
angular momentum i^. The three-term recurrence relation Eq. ()(M3|) has two inde- 
pendent solutions, which respectively behave for large n as 

An) (const.) (-g)" 

~ r(n + (a + /3 + 3)/2-s) ' ^ ^ ^ 

^[2) ~ {const.)Cr{n + {a + + l)/2 + s). (C-18) 

The first one, ^("j*) is the minimal solution, and the second one, ^(2)' ^ domi- 
nant solution, since lim ^^??M^o^^ = 0. In the case of the dominant solution these 

n— »oo y'^) 



coefficients A^^} increase with n, and the series ()(Mn|l diverges for all values of 



'(2) 



X. 



In the case of the minimal solution this series converges. Thus, we have to choose 
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in the series expansion ()(M0|) . For a general sEgmi^,), ^(\) does not satsify 
the relation HC-12|) . Hence, the requirement to satisfy this condition determines the 
descrete eigen values gEimiC)- 

As a practical way to obtain A^^^^ as well as sEim{C), we introduce 

Rn = Ln = -—Y- (C-19) 

^(1) ^(1) 

The ratio i?„ can be expressed as a continued fraction, 

Q,{")^("+1) a^n+l)^{n+2) 

We can also express in a similar way as 



Li; 



^(n) c,(«-i)7W a(«-2)7(-i) 0(1)7(2) 0(0)7(1) 



^(n)_ /3(n-l)_ /3(n-2)_ ^{1)_ ^(0) 



(C-21) 



This expressions for i?„ and L„ are valid if the continued fraction ()C-20() converge. 
(Notice that the last step of Eq. HC-2ip is not a continued fraction, but just a rational 
function.) By using the properties of the three-term recurrence relations, it is proved 
that the continued fraction ()C-20|) converges as long as the eigenvalue sEtmiO is 
finite. 

Dividing Eq. (RTTal) by the expansion coefficients s^^^, we obtain 

/?(") + a(")i?„+i + 7(")L„.i = 0. (C-22) 

We replace Rn+i and Ln-i by Eqs. (1^^201) and (IH^ . Then we can determine 
the eigenvalue s-E'^m as a root of Eq. (|C-22p . There are many roots, and the above 
equations for all value of n are equivalent. In practice, however, we truncate the 
continued fractions at finite lengths. In this case the most efficient way is to choose 
the equation with n = n£ := e-{a + P)/2. With this choice all terms in Eq. (l(T22l) 
become 0{^'^), and the length of the continued fractions that we must keep to achieve 
a given accuracy goal is the shortest. 

As was done in Fujita and Tagoshi's paper, in general, we can adopt Brent's 
algorithm^*^) in order to determine sEim{Cj- However, when |^| is not large, we can 
derive an analytic expression for gEimiO- The result is 

sEUO = ^(^ + 1) - ^^^e + [H{i + 1) - H{e) + 0{e), (C-23) 



with 



^(^^ - (2£- 1^3(2^ + 1) • (C-24) 
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After we obtain the eigenvalues sEimiOy can easily determine all the coef- 
ficients. The coefficient with n = is usually the largest term. The ratio of the 
other terms to the dominant term, i.e. A^^^/A^^^\ can be determined in the most 
efficient way with a minimal error due to truncation using Eqs. ()('-20|) and ()(y-21|) 
for < n < and n> ni, respectively. 

The coefficient of the leading term A^^^^ {sA^^') is determined by the normal- 
ization condition. Since ()C-4|) represents the same eigen function, the series HC-11|) 
should converge for the same eigen values sEimiOi constituting the minimal solution 
of the recurrence relation Eq. ()(M5|) . As in the case of {^l^-j*}, we have 



B 



in) 
(1) 



(n-1) 



B 



(1) 



B 



(n) 
(1) 



/3(")- 



a 



j3{n+l)_ p{n+2)_ 



a 



a 



(„_2)~(n-l) 



0(1)7(2) 



a 



(0)^(1) 



B 



(n+l) 
(1) 



I3(n)_ p{n-l)_ p{n-2)_ 



(C-25) 



(C-26) 



From these equations, we can determine the ratios of all coefficients, B^^^/B 

Now, we determine the values of the two coefficients A^^^^ and B^^'^ that deter- 
mines the overall normalization. Since Eq. ()C-7|) must hold for any value of x, we 
can set x = 1 in it to obtain 



^(1) • 



,B 



em ' 



r{n + a + 1) 



eM'^OsAt\OYl 



sA^;2iO nn + a + 1) 
A'-Z\Onn + l)na + l) 



n=0 s^im 



On the other hand, from the normalization condition (|C-2|) . we find 



dx 



1 



1 + x 



f3 00 



>(n2)p(a,/3)|. 



(C-27) 



l.(C-28) 



ni=0 n2=0 

Because the Jacobi polynomials are orthogonal, we have 



dx 



-1 



1 



X 



1 + x 



/3 



P(^'^)(x)p(^'^)(x) 



2r(n + Q + l)r(n + /3 + l),5„i, 



"2 



{2n + a + P + l)r(n + l)r(n + a + /? + 1) ' 
Then, Eq. (ITT^ reduces to 

2r(n + a + l)r(n + /3 + l) 



(C-29) 



E 

n=0 



An) 



.A 



(ne) 
tm 



,B 



(ne) 



(2n + a + /3 + l)r(n + l)r(n + a + /? + 1) 



Ant) Tjine)' 



(C-30) 
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Combining Eqs. ()(>27|) and ()('-3fl|) . we can determine the squares of and 



Bf'^\ Finally, we fix the signatures of s^)^' sB 



£m 

so that sS^^{x) reduces 



to the spin- weighted spherical harmonics in the limit ^ ^ 0. 
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